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I. INTRODUCTION 
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We present methods to compute maps of CMB fluctuations from high resolution cosmic string 
networks using a fuU Boltzmann code, on both large and small angular scales. The accuracy and 
efficency of these methods are discussed. 
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^ I The potential role of cosmic strings and other topological defects in cosmology has been the subject of considerable 
I interest for well over two decades (for a review see ref. [l]). Perhaps the most exciting prospect would be the detection 
' of their distinct observational signatures in the cosmic microwave sky. Cosmic strings, for example, are expected to 
create line-like discontinuities in the CMB temperature pattern, whereas other defects such as global monopoles or 
CN ■ textures create 'hot spots'. Their discovery would provide unprecedented information about the nature of unification 
^ " in the early universe, while their absence from the CMB would significantly strengthen constraints on a wide range 
of models. 

To many, the publication of the BOOMERanG results 0, in particular, signaled the demise of topological defects 
QQ . in cosmology. Indeed, the detection of an acoustic peak around £ ~ 200 was seen as evidence that primordial adiabatic 
perturbations were the seeds for large-scale structure formation, a view that has been strengthened with the apparent 
CsJ , resolution of further peaks (see also 0, 0, 01 )■ However, the presence of defects is not incompatible with inflation and 
• post-BOOMERanG analyses, such as refs. H,l3i concluded that they could not be ruled out. Current data allows 
' defects to play a significant (but subdominant) role in large-scale structure formation. In this sense, it is of great 
O i' importance to accurately characterize nonCaussian signals from strings, as they are likely to provide the only direct 
method of detection. 

In this paper we will detail the methods that we have developed to create full-sky and high resolution CMB 
maps generated by cosmic defects or any other 'causal' or 'active' sources. First, in f|n]we detail the large set of 
^ , perturbation equations that have to be solved, following this in illlll with a discussion of the treatment of the source 
terms which distinguish this analysis from that for inflationary fluctuations. In ill VI we then discuss efficient numerical 
implementation of CMB map-making using the analogue of Green's function techniques, without which the problem 
would not be tractable computationally. However, we complete this introduction by discussing previous work on 
?-j [ cosmic defects and the CMB, pointing out its relationship to this paper. 

Some of the earliest work featured analytic results obtained for simple string configurations 0, ^3 ■ Such exact 
solutions are important for testing computational methods. However, although these analytic results are interesting, 
numerical simulations are essential to obtain accurate quantitative predictions in more general contexts. The main 
drawback of numerical results in this context is their limited dynamic range, restricted by the light-crossing times 
of cosmic defect simulations. All-sky (large angle) CMB maps can be generated and have been used to obtain the 
normalization of the power spectrum to COBE, but their angular resolution has been poor. On the other hand, small 
angle maps permit the very important characterization of nonCaussian signals due to defects. Given the prospect of 
high resolution all-sky observations from the MAP and Planck satellites, ideally one would aim to compute all-sky 
defect maps of corresponding resolution, but computational resources remain insufficient for this task at present. 

Probably the earliest attempt at computing realistic CMB patterns generated by defects was that of Bouchet, 
Bennett & Stebbins j^. They employed a flat-space formalism to calculate the CMB temperature AT/T in the 
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direction n, solving the metric perturbation equations using Green's functions G(k, t, t') schematically as: 

-(n,k,t)cx / G^,(n,k,t,i')e^.(k,t')dt', (1) 



AT,. 



where 0^^, is the energy-momentum tensor of cosmic strings. Their methods neglected many effects, notably the 
presence of baryons and the expansion of the Universe, concentrating solely on the integrated Sachs- Wolfe effect. Pen, 
Spergel & Turok computed all-sky maps (at COBE resolution) produced by different global defects including an 
approximate treatment of CDM, baryons and radiation and their work was extended on intermediate angular scales in 
ref. 13J. COBE resolution maps generated by local cosmic strings were presented in using the AUen-Shellard (AS) 
string code j^^. The power spectrum for this map was evaluated for £ < 20 (using an ensemble of 192 realisations) 
and they inferred the string linear energy density to be G/z/c^ = I.OS^q ^q x 10~^. 

The most recent work on CMB fluctuations in the presence of causal seeds have made use of full Boltzmann codes 
(see section^, thus including all the relevant physics (to first order). The AS string code was employed again in 
the full Boltzmann analysis in ref. |l6j | , in which power spectra were computed from the brightness distribution, thus 
bypassing the maps. Power spectra were computed from simulations of different cosmological epochs and provided 
clear evidence of the importance of vector and tensor modes in these models, as well as the apparent absence of 
strongly defined acoustic peaks. 

An alternative line-of-sight approach was used in ref. ^-^i^ ^^^o later in refs. [TsLfigl l. to calculate power spectra 
for global defects. Here, the idea was to use unequal time correlators (UETCs) of the defect energy-momentum tensor 
(approximated by an expansion in eigenvectors) as the source for the perturbation power spectra. In principle, the 
method greatly extends the available dynamic range by exploiting the scalability of the correlators during defect 
evolution. However, while scalability is approached asymptotically in the radiation and matter eras, during the 
important radiation-matter transition the UETCs must still be calculated from large simulations bridging this time 
period. The line-of-sight method has also been used by (l^, El who employed an ensemble of toy model realisations 
of a string network and averaged the power spectra. The line-of-sight method can be used to compute maps as well: 
Simatos & Perivaropoulos |22 | modified it using a more general expansion of plane waves to accommodate for phase 
differences in a toy model for wiggly strings. However, while the method is phenomcnologically interesting it was 
necessary to make a number of assumptions about the string perturbation phases. 

It is important to note that none of these methods is perfect, and in some sense, they are complimentary: The direct 
approach developed further here, solving the full Boltzmann equations on a three-dimensional grid, provides reliable 
high resolution CMB maps. However, the UETC method with a greater dynamic range provides a more extensive 
view of the angular power spectrum. 

II. COSMOLOGICAL PERTURBATIONS 

In this section, we shall derive the equations that describe the evolution of first order perturbations in the metric 
and the energy-momentum tensor of matter fields in the presence of causal seeds such as cosmic strings. The formalism 
used is similar to that of [12| . except here, the full Boltzmann equation for relativistic matter (photons, neutrinos) 
is used. The treatment of the Boltzmann distributions presented here follows on from the approach used for scalar 
modes by '2^, but is extended to vector and tensor modes. The treatment of photon polarization follows the approach 
used by 2A\ for scalar and tensor modes, but is extended to vector modes. Similar, though not as general, methods 
were used by p^. where scalar and vector modes were treated in a gauge invariant formalism. 

A. Scalar - Vector - Tensor Decomposition 

In this paper, we will restrict ourselves to flat FRW models because present cosmic defect simulations are restricted 
to these. However, this is a significant simplification which enables us to expand all perturbations in terms of Fourier 
modes. In Fourier space, a tensor quantity Tij can be decomposed into scalar T, , vector and tensor components 
in the following way: 

T,,(k) = Is^.T + ihh - ^<5y)T^ -I- {kTj + kjTY) + i;^ . (2) 
This is analagous to the manner in which vector quantities can be decomposed into scalar and vector components: 

V,{]^) = hV^ + . (3) 



3 



In the above, vector and tensor components are tranverse that is, kiV^^ = kiTY = kikjT^j = and tensor components 
are, in addition, traceless T^i = 0. It is also useful to express vector components in an orthonormal basis ei, 62 with 
Gi X 62 = k, so that 

V = V^k + V^e^ + y/e2 . (4) 
We can also construct a basis for the tensor components out of ei and 62 by defining the following two matrices: 

M+ = ei (g) ei - 62 (8) 62 

Mx = ei (g) 62 + 62 (8) ei . ^ ' 

Tensor components can then be written 

rJ = Tj(M+),,+Tj(Afx)„-. (6) 

B. Einstein Equations 

We define perturbations of the conformally flat FRW metric as: 

gtiu = a^{'n){Vt^u + h^,,), (7) 

where rj^^ is the Minkowski metric and \h^j,v\ <C 1. We will work in the synchronous gauge in which /iqo = ^oi = 0- 
This gauge is not completely specified. This will result in extra "constraint equations" from the Eintein equations, to 
ensure that all degrees of freedom are specified. Perturbations in the energy-momentum tensor are given by: 

5n = -5p + Ql 

5Tf = {p + p)v, + Ql (8) 
ST} ^ 5p6] + pY^j + , 

where is the energy momentum tensor of the causal stiff source and contains the anisotropic stresses of 
relativistic matter. For the velocity terms, Vi and 0oi, we use slightly different variables to those defined in section lTl Al 

vY = , vY = ivY ■ 

These new variables enable us to write the evolution equation without explicit i's, which is easier to implement 
numerically, as will be discussed in section Hvl 

The equations for h^j^^, are then obtained by substituting equations ((TJ and (jSJ into the Einstein equations, using 
the homogeneous part for the background metric. The 00 and ij components of the equation 

5{Rf,, + Ag^,) = 87rG5(T^, - ]^g^,T^) (10) 

then lead, to first order, after transforming to Fourier space and decomposing into SVT components, to the following 
equations of motion for h : 



a e 

and the 00 and Oi components of the equation 

1 

2 

lead to the following constraint equations for h: 



h+^h^ -8TTG{a^{Sp + 3Sp) + Goo + O) , 
'hS 2f/i^ + ^h- = IGTrGia^pE^ + 6^) , 

hf + 2^hY = mnGia^pJ^Y + ^Y) . 
hT + 2^h^ + k^h^ = 16^G(a2pS^ + Q^) , 



(11) 



S{R^. - -Rg^. - Ag^,,) = 87:GST^, (12) 



k^h" + 3f /i = 2ATrG{a^Sp + 600) , 

=24TTG{QD-a^{p + p)0), (13) 
khy = 16ttG{VY - a'{p + p)dY) > 
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where i = 1, 2 and e = +, x and wc have defined h~ = h — . 

Energy conservation T'^''.^ = leads to the foUowing equations of motion for matter and radiation perturbations: 

S = -{l + w){9 + lh)-3^{cl-w)S, 
= - icl)e + ^kH + liq^fc^E^ , (14) 

^r = -f(i-3a«r-T^^fcsr, 

where w = p/p, = 6p/6p is the sound speed squared and 6 = 5p/p. These equations are vahd for uncoupled fluids. 
And covariant energy conservation with respect to the background G||^ leads to 

eoo = -f(eoo + e) + ez) 
K = -K-PY + kef • 

In what follows, we shall drop the tilde on the new variables to render the equations more legible 

C. Relativistic Matter 

To treat photon perturbations, wc derive the equations of motion for the Stokes parameters, which describe polarized 
light: the intensity /, the orientation of the polarization ellipse Q and U and the ratio of its principal axis V. It is 
convenient to work with the perturbations normalized to the average intensity Iq = p-y/^Tr: 

/ = /o(l + A/), 
P = loAp , 

where P stands for Q, U or V. We start with the general transfer equations for polarized light: 

A/ = +ikiiAi + 2hijhihj = f (Af - A/ + 4n • Vb) , 
Ap + ikiiAp = f{Af - Ap) , 

where /x = n • k, f = aaxng is the differential Thompson cross-section, Vb is the baryon velocity and the S denotes 
scattered quantities as measured in the comoving frame and the vector Ap has Ag, A^/ and Ay as components. 

^From the start, it is convenient to notice that V is always uncoupled to the other parameters as it cannot be generated 
through Thompson scattering. So we can set it to zero without loss of generality. We shall split the perturbations 
into scalar, vector and tensor parts: Ax = A|- + A^ + Aj, where X = I, Q or U. The treatment for massless 
neutrinos is identical to that for the intensity of photons, except that all terms involving f are zero. 

1. Scalar perturbations 

For scalar modes, U is also uncoupled from I and Q, and hence can be set to zero. The contribution from the 
metric is 

n,n,Aff = ^ + (m' - l)h' 

and the "polarization term" is given by 

where the relevant block of the "scattering matrix" M^{p,,fi') is given by 

/ 3 - - ^2 ^ 3Ai V" 1 - - 3/i2 + 3M2/i'2 \ 

V 1 - 3/x'2 - + 3^2^/2 3 _ 3^/2 _ 3^2 ^ 3^2^,2 j 



(15) 



(16) 



(17) 



(18) 
(19) 

(20) 
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To integrate out the angle dependence, we expand the perturbations in Legendre polynomials: 

= E(-^)'(2^ + ^)^xiPe{n) ■ (21) 



The equations then read 



= f(A5^ + W-Af-iP2(/im, (22) 
A§ + ifcM = f(-A^ + (1 - iP2(M))n^) , 



where 11^ = A?, + A^n + AS. 



2. Vector Perturbations 



Unlike scalar perturbations, U does not decouple from I and Q in the vector case, but this can be dealt with very 
easily. First we consider the contribution from the metric: 

hihjh'i^'*°'' = 2n^/l^{hY cos ip + hX sin ip) , (23) 
where cosy = n • ei and sin^? = n • 62- The polarization term is given by 



A^^(/x) = ^ / yr^v'T^M^CM, p', ^)AUf^')dn' , (24) 



where "& = (f' — ip, and the scattering matrix M^(/i, 1?) is given by: 

(2/x^'cosi? 2^^' cosd ^siwd 
2/^/i'cosi? /isinz? . (25) 
— 2/x'sin-(^ — 2/x'sin-!? cosz9 

Unlike the scalar case, there is a dependence on the azimuthal angle. To eliminate it, we introduce new variables 
defined as follows: 

Ay = -iy /T^ {AY^ cos if + Ay2 sin if) 

A^ = fl ^/l^ {A^^ cos if + A^2 (26) 
= VT^^(- Aj;i sin (f + A^2 cos cp) . 

With these variables, the equations for each components decouple and depend only on n so that we can decompose 
the new variables into Legendre polynomials. The equations then read 

Aj' + iknAj' + 4ifihY ~f{AY' - AvYi + ifill^') , 

A^^ + ikiiAl^ = -f(A^» - , (27) 
A§' + A^« = 0, 

wnere ii - jq'^h + io^/3 + M^QO 4'^Q2 + 35^Q4- 

3. Tensor Perturbations 

Again, the Stokes parameter U cannot be set to zero, but, as with vector perturbations, we will show that can be 
treated with Q. The metric term is 

fnfijhfj''"''' = {1 - n'^Xhl cos 2(p + hi sin 2y) . (28) 

The polarization term is given by 

a^^(m) = ^J m^(m, m', , (29) 
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where the scattering matrix M'^ fj,' ,-d) is given by 

/ S_(Ai)5-{/i')cos2i9 -S-{^l)S+{^l')cos2■& S-ifi) sin 2^ \ 

- S+ (ti)S -ill') COS 2-0 5+(/i)5+(At')cos2i? ii' S+i/i) sm2^ , (30) 
\ 2fiS-{fj.')sin2-d -2fiS+{fj.') sin2-d 2 fifi' cos 2-d J 

where S± (/i) = 1 ± . 

As with the vector case, we can ehminate the ip dependence by making a change of variables. In the tensor case, 
these new variables are defined as follows: 

= (1 - /i2)(Aj+ cos2(p + Aj"" sm2ip) , 
Ag = {1 + n^){A]^+ cos2ip + A^"" sin 295), (31) 
A^ = 2^(- A^+ sin 2^ + A^"" cos 2ip) . 

Exactly as for vectors, the equations for each polarization decouple and depend only on fi so that we can expand 
them in Legendre polynomials and carry out the integration. The resulting equations are: 

Ap + ifc/iAp + 2/if = -t(aP + n^^) , 

A5^ + ifc/.A5^ = -f(A5^-n^'), (32) 



aS^ + ap = o, 



whpTP n"^^ - J-A'^'= 4- iA^'= 4- -2-A^« - 2a^'= -I- SA^« - -2-A'^'S and ^ - + or x 
wnere li - jg^i^g + ^^^2 + 70^/4 s^QO + 7^Q2 7o^Q4 ^^^^'^ - + or x . 



^. Physical significance of moments 



If we use the decomposition (|21ll and the recursion relation 



2^+1 

we obtain the following equations for the moments: 



/^^^(m) = 7T7^m-l(^) + ii+l)Pi+l{fi)) , (33) 



Af„ = -A:Af,-|A, 

Afi = -|(2Af2 - Afo) - r(Af, - A^,) , 

Af2 = |(2Afi - 3Af3) - f(Af2 - ^,U^) - ^.h'^ , (34) 
Af, = 2FPT(^Af,_i -{£+ l)Af,^ J - fAf, , £ > 2 , 

^Qi - 2JTt(^A3^-i - + l)AS,+i) + r(-A§, + n^(<5,o + i^fe)) , 



AS = -fcA]r/-f(Aro^ + 40, 

Ar^* = liAY,^ - 2AY^) - Iky - r{AY,^ + n^o , 

Ar/ = 2?tt(^^«-i - + i)^)'^Vi) - ^Ar/ , ^ > 1 , 

A^l = ^A^^U-i l)A^^+i) - rA^^ + iU^%e , 



Ap = -A:AP-2/.J'-f(Ap + n^^), 

AP = 5^(Mf/_i - (£ + l)Aj;VJ - rAjl , £>0, (36) 
Agl = - + l)AS^,+i) - tA^^, - rn^^6,o • 

These moments are related to the energy momentum tensor of photons and (massless) neutrinos. First, let us 
consider the energy density. We have that 

SToo^-6p^-£j Ajdn, (37) 

so that 



^10 ■ 



(38) 



7 



Secondly, the velocity perturbations 

r)7n,. = ( n -i- r)"l?),- — 

An 



STo^ = (p + p)v^ = ^ I Ain,dn , (39) 



which yields 



Finally we have the anisotropic stresses 



which yields 



E5 = -3Af2, 



sr =|(A,7 + A)'3»), (42) 

— 5^/0 + 7^/2 + 35^/4 ■ 

We now have an open hierarchy of equations for all the components of the energy-momentum tensor. Numerically, 
it has been found that truncating the expansion by simply setting the moment to propagates a sizeable error 
back to the first moment. It has also been found 23] that the following expression gives excellent results: 



A,_^ = A:A,„„_i - ^-1^^ + fj A,_^ , (43) 

is the moment at which the series is truncated. The above expression is good for all photon and neutrino 
equation hierarchies. 



D. Non-Relativistic Matter 



For non-relativistic matter, it is not necessary to consider a full phase space expansion, as all but the first few 
moments are totally negligible. Our starting point will be the equations derived from the conservation of the energy- 
momentum tensor p4|l . 

CDM does not couple with other types of matter (except gravitationally), so we can immediately use the conservation 
equations. Here w = = 0. Also, in the (comoving) synchronous gauge, CDM has Vc = = 0, so H14|l simplifies to 

'5c = -^. (44) 
For baryons, these conservation equations H14|l become 



+ 1) 

^hi — a "hi 7 



' -^P6b (45) 



where we have used the fact that w, cj. ^ 1. However, baryons interact with photons through Thompson scattering. 
So we must correct the last two equations for momentum exchange between the two fluids. The equations for photons 
obtained in section ITl C 41 read 



.^-fc2(i<5^-il]S)+f(0f,-0^) 

By comparing with the conservation equations with w ^ — ^, we see that there is an extra drag term f{vi, — v^). 
Hence, for momentum to be conserved, we have to add the term |^T(w-y — Vh) to the previous two equations, which 
then become 

(47) 



where we have defined R = ^f-^ . 
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E. The Sachs- Wolfe Formula 



The temperature anisotropy formula is derived by considering perturbations in the photon energy along the 
unperturbed path X'^ = h'^ri. Here n is the line of sight direction (i.e. pointing in the direction opposite to the 
photon's travel). For a proof, see e.g. |l2ll25l |. In the synchronous gauge, the temperature fluctuations are given by: 

— = — -v^-n- - / h.jn.Ujdr] . (48) 

The above formula was derived under the assumption of instantaneous recombination. To treat the finiteness of 
the surface of last-scattering, we integrate the expression for the temperature fluctuations over the probability of 
free-streeming for a photon, e~'^dT: 

— =y —re dr] = J ^re [jj ^ J ^ 2^ h.jn.rij j dr] , (49) 

where the ISW term was obtained by integrating by parts and setting the surface term to zero because the visibility 
function is very sharply defined around the time of recombination and hence is utterly negligible a,t rj = or today. 

The calculation of the thermal history of the universe for an arbitrary set of cosmological parameters was achieved 
with an integrated package which will be described in more detail elsewhere p^ . The Friedmann equations were 
solved simultaneously with the ionization rate equations for hydrogen and helium. The results were compared for 
accuracy against RECFAST '27^ for which the code provided an independent check '26^. These computations start 
very deep in the radiation era and end today in order to create high accuracy tables from which the relevant quatities 
for the Boltzmann evolution, such as the opacity and visibilty function, are later interpolated using a cubic spline. 



III. COSMIC DEFECT SOURCE TERMS 



The perturbation source terms in the Boltzmann evolution are given by the Fourier transform of the cosmic defect 
energy-momentum tensor, decomposed into scalar, vector and tensor components. The code can accept the energy- 
momentum tensor of any set of 'active' sources with appropriate initial conditions, whether these are cosmic strings, 
global defects or other more exotic phenomena. We have experimented with inputting from global defect simulations, 
but the focus of our attention here is on local cosmic string simulations which inherently can achieve far higher 
resolution and greater dynamic range. 

The cosmic string simulations were performed using the AS code [isj . for which the methods employed and key 
results have been described in detail elsewhere. The strings are approximated by a two-dimensional worldsheet defined 
by x'^{a, t) — {t, Xs(cr, i)), where the position is a function of the two coordinates a (spacelike) and conformal time 
t in an unperturbed FRW background. We can impose the condition that the velocity x is transverse to the tangent 
vector x^ = dx^/da along the string, that is, x^ • x'^, — 0. The strings are evolved by splitting the equations of motion 
for the strings into their characteristic left- and right-moving modes, which are damped slightly by the expansion 
of the universe. The key points to note about the string network simulations are that above a minimum resolution, 
energy conservation is accurately satisfied during the numerical evolution to within a fraction of one percent, and over 
a dynamic range approaching an order of magnitude in conformal time (that is, several decades in redshift). 

The energy-momentum tensor of the strings is given by 

e^'-V^^/.,! da{ex^x:-e-'xTx'nS^'\x'^-K:{a)), (50) 

1/2 

where /i^ is the linear energy density of the string, e = (x'^/ (1 — x^)) , and g is the determinant of the background 
FRW metric {\/—g = a^). All the components of 0^i/(x, rj) were calculated at each point on the string network and 
then interpolated onto a high resolution grid. Since there is an implicit differentiation of the spatial distribution of (|5U|) 
in the scalar- vector-tensor decomposition described in section FlI Al the interpolation must be sufficiently smooth for 
the decomposition to remain well behaved. Simple cloud-in-cell interpolation, for example, produced poorly controlled 
results, so we chose a two step approach with a higher order scheme. The first step was to use the triangular shaped 
cloud (TSC) interpolation involving the 27 nearest neighbour points of the string segment. The weight function is 
given by 



FIG. 1; Positive and negative isosurfaces of scalar, vector and tensor components in real space created from the decomposition 
of the energy- momentum tensor of a straight string. Clockwise from top-left Ooo (positive), Os (positive centre and negative 
sidelobes), (negative centre and positive sidelobes) and Qy (alternating negative and positive). Also shown are contours 
in a plane transverse to the string. 

where |a;| is the distance between grid point and string in the x-direction. This weight function is multiplied by the 
appropriate weights in y- and z-directions. The second step is to smooth the result further by employing a gaussian 
window function in Fourier space, that is, proportional to exp(— fc^/A;^) with fc^^ corresponding to about three grid 
points; the discrete Fourier transform employed the NAG routine C06PXF. Energy-momentum conservation was closely 
monitored and maintained to within a fraction of one percent by this process. 

The smoothing from the interpolation and filtering leads to an apparent loss of resolution, but this need not be 
important. A well-behaved decomposition can be obtained at very high resolution on a larger grid and then used 
at lower resolution for the Green's function integration. We note also that, in Fourier space, the decomposed scalar, 
vector and tensor parts of H5U|) calculated as in section III Al can be stored on only half the complex grid, because the 
Fourier transform of a real quantity satisfies /(k) = /*(— k), it is necessary to treat only just over half of the grid, 
that is, ix, iy — 0, 1, — 1 and = 0, 1, ...,N/2 (it is in fact possible to use exactly half, but the relationship 
between components is more complicated |28j'). 

Figure n illustrates the decomposed scalar, vector and tensor parts of (|50|l for a straight cosmic string calculated 
as in section Hi Al The scalar Oqo is localised in real space, but the scalar and vector projection operators acting on 
the tensor components Qij yields the apparently non-local results shown. This demonstrates that the decomposed 
components must be calculated and evolved with great care; the final temperature pattern must reflect the locality 
(or causality) of the sources which generated them. 

IV. NUMERICAL IMPLEMENTATION 

The various evolution equations derived in section^ can be written as three systems (one for scalar, vector and 
tensor) of the following form: 



— = A{k, f])y + q{k, 7]) , 



(51) 
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where y and q are vectors of dimension n„ar, the number of equations to solve and ^ is a real square matrix of the same 
order. The components of y are the metric perturbation, CDM and baryon density contrasts and peculiar velocities 
and photon and neutrino moments and those of q are the defect source terms. Given that there is a continuous 
contribution at all times from these sources throughout space, we cannot simply solve this problem by projecting our 
initial conditions forward to the present time with transfer functions (as in inflation). However, rather than directly 
solving the differential equations at every grid point, it is much more efficient to employ a Green's function approach. 

In order to take advantage of the fact that the evolution equations depend only on the magnitude of the wavevector 
and not its orientation, we proceed to solve ()51|l by constructing the fundamental matrix Y of the system (e.g, see ^29^ ) 
which satisfies 

—Y-AY 

dr, ' (52) 

y(0) = 1 , 

where 1 is the identity matrix. Now, given some initial conditions 

y(0)=c, (53) 

the solution to H51() is 

y{v) = y{v) (c + £ Y-\r,')<i{r,')dr/^ . (54) 

The above equation reveals a similarity between the matrix Y~^ and a Green's function. However, the latter nomen- 
clature is not strictly appropriate in this case, as one of the basic properties of a Green's functions is undefined in a 
first order problem. 

A. Constructing the Matrices 

Since the system ()51|l is first order, we need to reexpress the metric equations in first order form (all the matter 
equations are already first order). For the scalar metric equations, we use the two equations of motion, with the first 
constraint equation to replace the k^h~ term, to obtain the following first order equations for h and : 



--h~STTGa^{5p + 2,6p)-STTG{Qa^ + Q), 

(55) 



drj 

flh'S A, 

^^-{h~ 2hS) + 8^Ga2(2pI]^ - 6p) + 8^G(2e^ - Ooo) • 
dr] a 



The vector metric equation of motion is already in a form that is first order for h 

dh^ d 



+ WiTGa^pT.^ + WnCe'' . (56) 



drj a 

For the tensor modes, we use the standard order reduction, by considering h"'" and hT' as separate variables: 

dh^ 



(57) 

, = -2-h^ - k^h^ + IGTrGa^pS'^ + IGttGQ'^ . 
drj a 



drj 

dh^ 



To construct Y, we numerically solve (|52|l . which involves solving 3 times (scalar - vector - tensor) n^ar systems 
of rivar coupled equations for a chosen number M of the wavenumbers k (typically M ^ N for an grid). This 
contrasts with solving 10 systems (real and imaginary parts of 1 scalar, 2 vectors and 2 tensors) of n^ar equations 
for /2 wavevectors k if we were to solve H51f) directly at every grid point. Hence this Green's function approach 
represents a reduction of computing time by a factor of 5N'^ /Sn^ar- Since n^ar — 40, this factor is more than 1,000 
for a 256'^ box and 10,000 for a 1024'^ box. Of course, there are extra computations involved: the inversion of the 
matrices (for which we use the NAG F07ADF and F07ADJ Fortran routines) and the actual integration of H54|l . but 
this is insignificant compared with the time taken to solve the extra equations. In fact, this method is so efficient it 
is constrained by the amount of memory (or disk space) required. 
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FIG. 2: Time evolution of elements of the fundamental matrices for rj — Srjdec to Arid^c- These particular diagonal elements 
are those associated with ISW effects, i.e. the time derivative of the metric terms: scalar trace h (solid), anisotropic scalar 
(dot-dashed), which coincides with the vector and tensor ft^ (dashed). Also illustrated is the CDM density contrast Sc 
(dotted). The conformal time along the horizontal axis is given in Mpc 



To numerically solve the perturbation equations, we use dverk.f, a sixth order Runge-Kutta Fortran routine. 
Figure 12 shows the time evolution for some diagonal elements of the fundamental matrices. Our Boltzmann code was 
extensively tested for inflationary scenarios and found to be in excellent agreement with publicly available codes [2^ . 
The matrix elements are computed for fc = to A: = 1.8NAk/2 with a spacing of 0.9Afc, where Afc is the simulation 
grid spacing. The value of 1.8 was chosen to be sightly higher than the maximum possible value VS- This is to ensure 
that there is no need to extrapolate the matrix elements. This fc-spacing is sufRciently dense for the elements of the 
matrices and their inverses to be linearly interpolated for the appropriate value of k at every simulation grid point 
(see FIG. El). 

Since only a few components of the vector q are non-zero (two for scalars, one for vectors and one for tensors), 
only the corresponding rows of need to be stored. Also, since only 10 quantities are needed to compute the 
SW integral {S^, 9j, w^, h, , hX and /i"^), only the corresponding columns of Y need to be stored in principle. In 
practise however, all components y need to be kept to allow the computations to be made in several stages. There 
are two reasons why we want to do this. The main reason is that, even though all the components of the fundamental 
matrices are well behaved at all times, the evolution results in different components having very large ratios. This 
results in Y become non-invertible to machine precision. This happens more often at early times and for scalar modes. 
The solution H54|) is then computed by 

y(^) = Y^^[^) (yin,) + /;;+^ y(^) -Hv'Mv')dv') 

rji <r] < , i = 0, 1, ...Ustage 

y{m=o) = c, 

where the corresponding matrices y'*^ and 1"^*^ ~^ are computed from rji to ?7i+i. The second reason for wanting to 
allow the possibility of running the computations in several stages is more practical: this enables us to checkpoint the 
code in order minimize the effects of a system crash. 
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FIG. 3: Diagonal scalar matrix elements Yn, Yjj^^ and Y33 from the late matter era corresponding respectively to h (continuous 
and dashed lines) and Sc (dot-dashed). The normalisation is arbitrary and the k scale is 1.8 X 10"^ Mpc'^. 

B. Initial Conditions 

As our simulations start at a time much later than that of the phase transition that created the defects, we must 
specify initial conditions that are consistent with a fully formed network. In any case, the local string simulations 
used are based on an effective action and hence cannot be used to simulate the phase transition itself. 

One of the nice things about the solution (|54|l is that it is expressed as a sum of a term depending on the initial 
conditions and a second one depending on the defect sources. This provides us with an easy way to assess the relative 
importance of the initial conditions and hence to see if our results are sensitive or not to them. 

One possible choice is to set all gradient terms and time derivatives to zero in the constraint equations. This leads 
to the following initial conditions: 







= js. 


Ob 


— 9^ ^ 




r 







Boo 
a'^ip+p) 



a^p+p) ' 



(59) 



the last line being valid for all three components. This choice is consistent with setting the pseudo-energy tqo and 
pseudo-momenta rpi to zero, a set of initial conditions consistent with matching an instantaneous defect-forming phase 
transition to a homogeneous initial state (see, for example, [sol l8l|). 



C. Integration over the String Energy-Momentum Tensor 



The matrix elements are integrated over the energy-momentum tensor using the trapezoidal rule. The equations 
being real and scalar, the same matrices are used for both the real and imaginary parts as well as the two vector 
and two tensor components, so that the integration is repeated twice for scalars, and four times for both vectors and 
tensors. 

The grids of the quantities needed to compute the SW integral, 6-y , and hij are then completed at each timestep 
by complex conjugation and Fourier transformed back to real space. They are then projected onto the end of the 



13 



vectors ni{r] — rjo), where the fii are the pixel directions, using the inverse cloud-in-cell scheme: 

f{x, y, z) ^ Wx{wy{wJ{xi,yi, zi) + (1 - ■w^)f{xi,yi,Z2)) 

+ (1 - Wy){W:,f{xi,y2,Zi) + (1 - W,_)f{xi,y2,Z2))) 
+ il-Wx){Wy{Wzfix2,yi,Zi) + {l-W^)f{x2,yi,Z2)) 
+ (1 - Wy){w^f{x2,y2,Zl) + (1 - W^)fix2,y2,Z2))) , 

where the xi, X2, ■■■ define the grid element and the weights are given hy Wx = 1 — xi and similarly for the other two. 
The temperature fluctuations are then computed using (I49|l . This integration is also performed using the trapezoidal 
rule as the integrand varies very smoothly along each line of sight. 



V. TEST OF THE NUMERICS 



In order to test the validity of our methods, we compare the result of simulations with analytic results. One of 
the difficult aspects of such comparisons is that analytic results are mostly obtained for defect configurations (i) in 
Minkowski space and (ii) in the no matter approximation. In an expanding Universe, with a complete treatment of 
matter perturbations, these approximations would be relatively good only in the late matter era and on scales much 
smaller than the horizon. 



A. Kaiser-Stebbins Effect 



To illustrate this, we consider an infinite straight string moving in a direction perpendicular to the linc-of-sight. 
According to 0, this will produce a discontinuity in the temperature fluctuations. This calculation is done in the limit 
that a plane wave of CMB radiation is propagating in the direction of the observer. Objects behind the string receive a 
boost towards the plane in which the string is moving because of the gravitational effect of its deficit angle A — SirGfi. 
This means that the CMB photons that were behind the string when they cross its plane will be blueshifted, so that 
they will be hotter than those that were in front of the string. The magnitude of the discontinuity is : 

AT 

— ^SnGpivj, (61) 

where v is the string velocity and 7 is the Lorentz factor. 

This provides an excellent test of our map-making pipeline, particularly the vector modes, which cannot be compared 
with inflationnary results such as is done for the scalars and tensors in chapter '26'|. However, in a realistic expanding 
Universe, such as those studied in this dissertation, this effect is more complicated due to the presence of matter, 
the growth (or decay) of SVT components, the curvature of the microwave sky and the decelaration of the string. 
In addition, this idealised solution must be studied in a periodic box with causal effects due to the limited dynamic 
range. 

In the realisation shown in figure ^ we minimised these effects by considering photon propagation in the late 
Universe with no cosmological constant, starting at i] = 0. 34770 and ending at r/ = O.4I770, on a square patch of sky 
of 3.2° X 3.2° through a box of comoving size L — 0.069?yo- The string was initially moving at a velocity v = 0.9, 
which had redshifted to v = 0.825 as the photons crossed the plane of the string. These test simulations included 
no compensation, i.e. everything was set to zero initially. The temperature jump measured from the average of the 
plateaus behind and in front of the string is approximately 

AT 

— c 5SG^/c' , (62) 

whereas the theoretical fiat-space, vacuum (no matter) prediction for the given velocities should lie in the range: 

AT 

— ^(37- 52)Gm/c^ (63) 

This is in good agreement with the expanding Universe simulation, despite the additional physical effects involved. It 
appears that the effect of matter on the scalar and vector modes slightly increases the anisotropics over the dominant 
flat-space ISW effect. 
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FIG. 4: Temperature discontinuity induced by a straight moving string with an initial velocity v = 0.9: total (top-left), scalar 
(top-right), vector (bottom-left) and tensor (bottom-right). The colour scale is given in units of A = SttG^i/c^ 

B. Loop Anisotropies 

In this subsection, we show the results of small angular scale simulations of string loops of radius R = 0.51]. The 
size of the horizon at the beginning of the simulation is 0.5L where L is the size of the box. The simulations end 
at 77 « lOijdec- Hence, we are dealing with large loops smaller than the horizon, but initially perturbed by Hubble 
damping. 

The general loop solution in flat-space can be expanded in Fourier modes, with the gauge conditions yielding 
constraints on the relationship between the different modes. These can be solved if only the first few harmonics are 
considered. The Kibble- Turok solution |32j| for a loop of length L involves the first and third harmonics: 

x(C, t) — {ei [(1 ~ k) sin ct- -|- Ik sin3(T_ + sin cr+] 

—62 [(1 — k) cos cr_ + ift cos 3(7_ -|- cos cos cr+] (64) 
—63 [2k-^/^(1 — k)^/^ cos cr_ -I- sin(pcos(T+] } , 
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FIG. 5: CMB fluctuations seeded by a single loop: total (top-left), scalar (top-right), vector (bottom- left) and tensor (bottom- 
right). 

where a± — (2'!t/L)(^± are the left/right moving modes and the e.^ are the cartesian unit vectors. 

Figure |5l shows the CMB fluctuations caused by a perfectly circular loop {(p — k — 0). The loop has an initial 
radius just below the angular scale of the patch shown in figure |S1 The loop begins stationary and then accelerates 
up to a relativistic velocity as it collapses. At this point it is crossed by the plane of CMB photons where a clear 
Kaiser-Stebbins discontinuity is evident; note that the amplitude of the temperature increases from large to small 
radius as the velocity increases. There is a cancellation of the non-local scalar, vector and tensor modes in the interior. 
By the time of photon crossing, causal effects related to the velocity acceleration have not had time to propagate to 
the centre. During one period, the average velocity of a loop is < w >= l/\/2. 

Next we consider a loop with tp — tt/3 and k — 1/2, which does not possess the simple symmetries of the circular 
case. Its initial state is one of maximum radius and minimum velocity from which it then shrinks, accelerates and 
begins to rotate and form cusps. FigureElillustrates the total anisotropy pattern viewed from three different directions 
into the box. In the top right corner, there are two things to note: when the photons first cross the loop, there is 
little effect as it is moving slowly; whereas there is a strong discontinuity later when the photons pass the front of the 
relativistically collapsing loop. The bottom left corner shows a side view of the loop, while the bottom right shows a 
much more complex anisotropy pattern as cusps form viewed from above. 



16 




FIG. 6: CMB fluctuations seeded by a single (Kibble- Turok) loop (top-left): front view (top-right), side view (bottom-left) and 
top view (bottom-right). 



VI. DISCUSSION 

We have presented the full set of sourced evolution equations with the Boltzmann hierarchy necessary to study the 
gravitational effects of cosmic defects on the CMB. We have focussed attention on using cosmic string simulations 
as the evolving causal source terms in these equations. We have developed efficient numerical techniques for the 
Green's function computation of high accuracy maps from these string network simulations. We have also presented 
numerical tests of the full pipeline. More quantitative results from extensive supercomputer simulations will be 
presented in |^ , featuring large-angle and small-angle maps respectively. 
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